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The  problem  of  the  computational  time  reversal  is  posed  as  the  inverse  problem  of  the  determination 
of  an  unknown  initial  condition  with  a  finite  support  in  a  hyperboilc  equation,  given  the  Cauchy  data  at 
the  lateral  surface.  A  stability  estimate  for  this  ill-posed  problem  implies  refocusing  of  the  time  reversed 
wave  field.  Two  such  two-dimensional  inverse  problems  are  solved  numerically  in  the  case  when  the 
domain  is  a  quadrant  and  the  Cauchy  data  are  given  at  finite  parts  of  coordinate  axis.  The  previously 
obtained  Lipschitz  stability  estimate  (if  proven)  rigorously  explains  and  numerical  results  confirm  the 
experimentally  observed  phenomenon  of  refocusing  of  time  reversed  wave  fields. 


1  Introduction 


1.1  Statement  of  the  inverse  problem 


This  is  the  first  publication  in  which  the  problem  of  computational  time  reversal  is  solved  numerically  via 
the  solution  of  an  inverse  problem  for  a  hyperbolic  equation  with  the  Cauchy  data  at  a  lateral  surface. 
Consider  the  standard  Cauchy  problem  for  the  hyperbolic  equation 


utt  =  Lu  ( r,t )  £  Mn  x  (0,T), 


(1.1) 


u\t=0  =  <P,  Ut  |t=0  =  Vh  (1-2) 

where  L  is  the  elliptic  operator  of  the  second  order  in  Mn  with  coefficients  independent  on  t,  the  function 
ip  £  H2(Mn)  and  the  function  if  £  H1(Mn).  Hence,  the  Cauchy  problem  (1),  (2)  has  unique  solution 
u  £  H2(M.n  x  (0,T)).  The  initial  conditions  ip  and  if  are  assumed  to  have  a  finite  support  D  C  Mn, 


sup<£>  C  D,  sup  ip  C  D. 


(1.3) 


We  are  interested  in  the  following 

Inverse  Problem  1.  Let  T  C  Mn  be  a  hypersurface  and  Ty  =  T  x  (0,T).  Assume  that  one  of  initial 
conditions  ip  or  if  is  known  and  another  one  is  unknown.  Determine  that  unknown  initial  condition  assuming 
that  the  following  functions  h  and  g  are  given 


u 


rT 


h, 


du 

dn 


(1.4) 


Below  we  call  the  problem  of  recovering  of  the  function  ip  “the  (^-problem” ,  and  the  problem  of  recovering 
of  the  function  if  “the  t/l-problem” .  Since  (1.4)  is  the  Cauchy  data,  then  the  Inverse  Problem  1  is  a  particular 
case  of  the  so-called  Cauchy  problem  for  the  hyperbolic  equation  with  the  lateral  data.  Uniqueness  and 
stability  results  for  this  problem  are  obtained  via  Carleman  estimates  and  can  be  found  in,  e.g.,  books  of 
Klibanov  and  Timonov  [9]  and  Lavrentiev,  Romanov  and  Shishatskii  [12],  as  well  as  in  papers  of  Klibanov  and 
Malinsky  [7],  Kazerni  and  Klibanov  [5]  and  Klibanov  [6].  The  Holder  stability  estimate  was  established  in  [12] 
and  the  stronger  Lipschitz  stability  estimate  for  bounded  domains  was  proven  in  [9,  7,  5].  Klibanov  [6]  has 
studied  both  ip-  and  i/i-problenrs  for  a  general  hyperbolic  equation  (1.1)  with  variable  coefficients  (including 
hyperbolic  inequalities)  for  n  =  2,3,  assuming  that  the  domain  D  is  either  a  quadrant  in  the  2-D  case  or  an 
octant  in  the  3-D  case.  In  the  case  of  the  quadrant  the  data  (1.4)  were  given  on  finite  parts  of  coordinate 
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axis,  and  they  were  given  on  finite  parts  of  coordinate  planes  in  the  case  of  the  octant.  Using  previous 
results  of  [9,  7,  5],  he  has  proven  the  Lipschitz  stability  estimate  for  this  problem,  has  shown  its  connection 
with  the  refocusing  of  time  reversed  wave  fields  and  has  proposed  a  convergent  numerical  method.  Note 
that  the  Lipschitz  stability  in  an  unbounded  domain  such  as  quadrant  or  octant  is  rather  surprising,  since, 
unlike  the  bounded  domain  case,  a  large  part  of  the  energy  never  reaches  the  surface  (curve  in  2-D)  where 
measurements  of  the  wave  field  and  its  normal  derivative  are  taken. 

The  authors  are  aware  only  about  two  publications  in  which  numerical  solutions  of  similar  problems  were 
presented.  The  first  one  is  of  Klibanov  and  Rakesh  [10],  in  which  the  method  of  quasi-reversibility  of  Lattes 
and  Lions  [11]  was  adapted  for  the  solution  of  the  Cauchy  problem  (1.1),  (1.4)  with  L  =  A  =  d2  +  <92  in  the 
square  with  the  lateral  Cauchy  data  at  the  boundary  of  this  square  (it  was  shown  in  the  recent  book  [9]  that 
the  quasi-reversibility  is  a  particular  case  of  the  Tikhonov  regularization  method,  and  convergence  rates 
were  established,  also,  see  [6]).  A  quite  good  robustness  of  this  method  was  demonstrated  computationally 
in  [10].  This,  observation  goes  along  well  with  computational  results  of  the  current  publication  and  might 
likely  be  atributed  to  the  existence  of  a  priori  Lipschitz  stability  estimate,  which  is  the  best  possible  one, 
also  see  Section  6.  The  second  publication  is  of  Kabanikhin,  Bektemesov  and  Nechaev  [4].  They  have 
considered  the  inverse  (^-problem  in  the  2-D  case  assuming  that  L  =  A  =  92  +  92, 

D  =  {(x,  y )  |  x  G  (0,  a),  y  G  (-6,  6)}, 

T  =  a  and  T  is  a  part  of  the  y— axis, 

T  =  {(ar,  y)  \  x  =  0,  y  G  {-b  -  T,  b  +  T)}. 

A  uniqueness  result  for  this  problem  was  established  and  numerical  simulations  were  conducted  in  this 
reference. 

The  present  paper  is  motivated  by  [6].  We  consider  the  numerical  solution  of  the  problem  of  [6]  in  the 
quadrant.  To  do  this,  we  specify  the  Inverse  Problem  1  as  follows. 

Inverse  Problem  2.  Let  the  function  u(x,y,t)  be  the  solution  of  the  standard  Cauchy  problem  (1.1), 
(1.2)  in  M2  x  (0,T)  with  L  =  A  =  02  +  <92.  Suppose  that  in  (1.3)  the  domain  D  is  the  rectangle  located  in 
the  first  quadrant, 

D  =  {(£,  y)  |  x  G  (0,  a),  y  G  (0,  b)}. 

Let  T  =  Ti  U  T2,  where  Ti  and  T2  are  parts  of  coordinate  axis  (see  Figure  1): 

Ti  =  {x  =  0,  y  e  (0,b  +  T)},  T2  =  {y  =  0,  x«E(0, a  +  T)}. 

Given  functions  h  and  g  in  (1.4),  determine  either  the  function  ip  assuming  that  ip  =  0  or  the  function  ip 
assuming  that  <p  =  0. 

Using  (1.3)  and  the  finite  speed  of  propagation  of  the  wave  field,  we  arrive  at 

^|9n  =  °,  t  G  (0,  T),  (1.5) 

where  the  domain  U  is 

=  {(x, y)  |  x  G  (-a  -  T,a  +  T),y  G  (-6  -T,b  +  T)}. 

For  7  =  1,2  let 

riT  =  Ti  x  (0,T),  hi  =  h\ViT,  gi  =  g\TiT. 

We  specify  conditions  (1.4)  as  follows 


4 

II 

Ux\v1T  ~  9ii 

(1.6) 

r2T  ^2, 

%lr2T  —  ^2’ 

(1.7) 

2 


Figure  1:  The  geometry  of  the  problem 


1.2  Connection  with  time  reversal 

Borcea,  Papanicolaou,  Tsogka  and  Berriman  [2]  were  the  first  ones  who  draw  the  attention  of  the  mathe¬ 
matical  community  to  the  issue  of  the  computational  time  reversal  (the  mathematical  model  of  [2]  is  quite 
different  from  ours).  Bardos  and  Fink  [1]  were  the  first  ones  who  have  noticed  the  direct  connection  between 
the  problem  of  the  computational  time  reversal  and  the  Inverse  Problem  1.  Later  this  was  also  observed  by 
Klibanov  and  Timonov  in  [9,  8].  The  direct  linkage  between  stability  results  for  the  Inverse  Problem  1  and 
refocusing  in  time  reversal  was  first  observed  in  [6] . 

The  set  up  for  the  time  reversal  experiment  is  well  described  in  the  review  paper  of  Fink  and  Prada  [3]. 
A  point  source  at  the  moment  of  time  t  =  0  causes  a  wave  pulse.  As  a  result,  waves  propagate  in  M3.  At 
a  surface  S,  both  the  wave  field  and  its  normal  derivative  are  recorded  by  many  transducers  for  the  time 
period  t  £  (0,T)  (see  p.  R3  in  [3]).  Next,  the  wave  field  is  “send  back”  from  those  transducers  by  the 
principle  “first  in-last  out” .  Actually,  therefore,  “sending  back”  the  wave  field  means  that  experimentalists 
solve  experimentally  the  problem  of  finding  the  solution  of  a  hyperbolic  equation,  given  the  Cauchy  data 
at  S,  i.  e.,  given  the  wave  field  and  its  normal  derivative  at  S.  However,  the  initial  condition  at  {t  =  0}  is 
unknown  in  this  case.  Furthermore,  the  determination  of  this  initial  condition  via  refocusing  of  the  time 
reversed  wave  field  is  actually  the  main  “goal”  of  the  time  reversal  experiment.  Consequently,  a  direct 
3-D  analog  of  the  Inverse  Problem  1  is  solved  experimentally.  In  doing  so,  an  experimental  device  actually 
“constructs”  the  function  v(X,t)  =  u(X,T  —  r),  r  £  (0,  T),  X  £  R3  where  u(X,t )  is  an  approximation  for 
the  function  u(X,t).  The  major  practically  interesting  phenomenon  observed  in  these  experiments  is  that 
the  time  reversed  wave  field  refocuses  at  the  location  of  the  original  source.  In  other  words,  it  replicates  the 
5— like  function,  which  represents  the  source  of  the  original  pulse. 

It  is  an  interesting  mathematical  problem  to  explain  the  refocusing  phenomenon  rigorously  and  to  confirm 
it  numerically.  The  phenomenon  of  refocusing  in  the  quadrant  is  demonstrated  numerically  in  this  paper. 
In  terms  of  the  function  v{X,  r),  refocusing  means  that  this  function  at  t  ~  T  is  close  to  the  5-function, 
which  represents  the  original  point  source.  Therefore,  to  rigorously  prove  refocusing,  one  should  establish  a 
stability  estimate  for  the  Inverse  Problem  1.  Indeed,  a  stability  estimate  guarantees  that  the  approximate 
solution  u(X,  T  —  t)  should  be  close  to  the  exact  one  u(X,  T  —  r),  as  long  as  the  measurement  error  in  the 
data  at  the  surface  S  is  sufficiently  small.  Hence,  the  function  v(X,  r)  =  u(X,  T  —  r)  «  u(X,  T  —  r)  should 
be  close  to  the  original  5-function  at  r  «  T.  Consequently,  the  Lipschitz  stability  estimate  of  the  paper  [6] 
in  the  quadrant  ensures  both  a  high  degree  of  refocusing  in  the  quadrant  and  a  good  stability  property  of 
the  numerical  solution  of  the  Inverse  Problem  2.  We  observe  the  latter  in  numerical  experiments  (Sections  5 
and  6).  We  also  mention  that  previous  Lipschitz  stability  estimates  of  [9,  7,  5]  were  obtained  for  the  case 
when  the  lateral  Cauchy  data  are  given  at  the  boundary  of  a  bounded  domain  and  solution  is  sought  for 
inside  this  domain.  Thus,  they  imply  a  high  degree  of  refocusing  of  time  reversed  wave  fields  propagating 
in  bounded  domains. 
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1.3  The  Cauchy  problem  with  the  lateral  data  is  not  equivalent 
with  the  Inverse  Problem  1 

Although  the  Inverse  Problem  1  is  a  particular  case  of  the  above  mentioned  Cauchy  problem  (1.1),  (1.4)  for 
the  hyperbolic  equation  with  the  Cauchy  data  at  a  lateral  surface,  but  these  problems  are  not  equivalent. 
Indeed,  consider,  for  example  the  1-D  case.  Let  the  function  w(x,  t )  be  the  solution  of  the  standard  Cauchy 
problem 

wtt  =  wxx ,  (x,t)  G  M  x  (0,T), 

rc(x,0)  =  a(x),  wt(x,  0)  =  /3(x), 

where  a  G  H2(K)  and  [3  G  i/1(M).  Suppose  that  the  following  two  functions  fi(t)  and  /2(f)  represent  the 
lateral  Cauchy  data  at  {x  =  0} 

w(0,t)  =  wx(0  ,t)  =  f2(t),  te(0,T). 

By  the  D’Alembert  formula, 

x+t 

.  .  a(x  —  t)  +  a(x  +  t)  1  f  ....  .. 

w(x,t)  =  - - j  /J(£)d£. 

x—t 


Let  a(x)  =  /3(x)  =  0  for  x  <  0.  Then  we  obtain 

A<*>  =  er  + 1 /*«*■*  = 

0 

These  two  equations  are  not  independent  ones,  since  the  second  equation  can  be  obtained  from  the  first  via 
the  differentiation.  Therefore,  because  of  the  presence  of  two  unknown  functions  a  and  (3,  the  latter  system 
has  infinitely  many  solutions  if  f[(t)  =  /2(f)  and  has  no  solutions  if  /( (t)  /  /2(f). 

The  fact  that  functions  fi(t)  and  /2(f)  are  actually  not  independent  from  each  other  has  the  following 
explanation.  Suppose  that  the  hypersurface  T  in  the  Inverse  Problem  1  is  the  boundary  of  a  certain  domain 
G  C  Mn  and  D  C  G  (G  can  be  both  bounded  or  unbounded,  in  the  case  of  the  above  example  G  =  {x  >  0}). 
One  can  uniquely  solve  the  boundary  value  problem  for  the  equation  (1.1)  with  the  zero  initial  condition  at 
{t  =  0}  and  with  the  Dirichlet  boundary  condition  u |rT  =  h  in  the  domain  (Mn  \  G)  x  (0,  T).  Solution  of 
this  problem  uniquely  determines  the  function  g  in  (1.4). 

Another  explanation  of  this  example  of  the  non- uniqueness  is  that  being  given  for  t  G  (0,T),  functions 
fi(t)  and  j\  (t)  determine  the  function  w(x,  t )  uniquely  only  in  characteristic  triangles  {x  >  0,  x  <  t  <  T  —  x} 
and  {.x  <  0,  —  x  <  t  <  x  +  T}  in  the  (x,  t)-plane,  and  the  axis  {t  =  0}  intersects  with  these  triangles  only  at 
the  origin.  On  the  other  hand,  assume,  for  example  that  the  function  a(x)  =  0.  Then  the  function  w(x,  t ) 
has  the  odd  extension  w  G  H2(R  x  (0,T))  in  {t  <  0},  so  as  functions  fi(t)  and  /2(f).  In  both  cases 
characteristic  triangles  become  {x  >  0,  x  —  T  <  t  <  — x  +  T}  and  {x  <  0,  — x  —  T  <  t.  <  x  +  T}.  Since  the 
interval  {t  =  0,  0  <  x  <  T}  is  the  median  of  the  first  triangle  and  the  interval  {t  =  0  —  T  <  x  <  0}  is  the 
median  of  the  second,  then  the  function  (3{x)  can  be  determined  uniquely  in  this  case.  In  the  case  /?(x)  =  0 
we  obtain  the  even  extension  and  repeat  these  arguments. 

Thus,  the  uniqueness  of  the  Cauchy  problem  with  the  lateral  data  does  not  necessarily  imply  uniqueness 
of  the  inverse  problem  of  the  determination  of  initial  conditions.  This  means  that  the  title  of  this  subsection 
is  true.  This  discussion  also  demonstrates  that  it  is  unlikely  that  both  initial  conditions  ip  and  ^  can  be 
determined  simultaneously  in  Inverse  Problem  1  and  2,  at  least  in  the  case  when  T  is  the  boundary  of  an 
unbounded  domain.  This  justifies  statements  of  Inverse  Problems  1  and  2.  Note,  however  that  this  issue  is 
a  subtle  one.  Indeed,  if,  for  example  P  is  the  boundary  of  a  bounded  domain  and  both  initial  conditions  p 
and  are  sought  for  in  this  domain,  then  results  of  [9,  7,  5]  imply  that  both  these  initial  conditions  can  be 
determined  uniquely  and  simultaneously. 
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2  Numerical  method  for  the  direct  problem 


To  solve  the  direct  problem  (1.1),  (1.2),  (1.5),  we  use  the  method  proposed  in  [4].  The  zero  boundary 
conditions  (1.5)  allows  us  to  represent  the  solution  in  the  form  of  the  Fourier  sinus-series 


u(x,y,t)  =  J2  Unm  (t)Xn(x)Ym(y), 


„  ,  .  .  im(x  +  a  +  T) 

Xn[x)  =  sin 


!„(„)=  sin  ™(9  + 6 +  T) 


2  (a  +  T)  ’  ““  2(6 +  T) 

Substituting  (2.1)  in  (1.1),  (1.2),  for  unm(t)  we  obtain  the  following  problem 


7 r2n2 


^"nm  ^ nm  ^nm  > 


+ 


ir2m2 


nm  4 (a  +  T)2  '  4(6 +  T)2’ 

Hran( 0)  =  -pam  •  Unnl(0)  =  Ipnm ■ 

where  <pnm,  ipnm  are  respectively  the  Fourier  coefficients  of  the  functions  ip(x,y)  and  i/j(x,y), 

1 


^Pnm  — 


V^nm  — 


(a  +  T)(6  +  T) 

1 


(a  +  T)(6  +  T) 
The  solution  to  the  problem  (2.3),  (2.4)  is 


<p(x,  y)Xn(x)Ym(y)  dx  dy, 
ip(x,  y)Xn(x)Ym(y)  dx  dy. 


—  prim  COS  (uJnrnt ) 


i>r 


UJn 


Sin  (c Onmt). 


Taking  into  account  (2.1),  we  obtain  the  following  formula  for  solution  to  the  direct  problem 
u(x,y,t)  =  V]  ( tprim  cos  (unmt)  +  sin  (unmt))xn(x)Ym(y). 

This  formula  can  be  written  in  the  operator’s  form 

u  =  Aq, 


(2.1) 

(2.2) 

(2.3) 

(2.4) 

(2.5) 

(2.6) 

(2.7) 

(2.8) 

(2.9) 


where  A  is  a  linear  operator,  q  is  the  vector  of  Fourier  coefficients  of  the  function  <p  in  the  99-problenr  and 
of  the  function  7/1  in  the  ^-problem. 


3  Numerical  method  for  the  inverse  problem 

Let  /  =  (h,g)  be  the  vector  function  representing  the  data  (1.4)  for  the  inverse  problem.  The  Inverse 
Problem  2  is  a  linear  inverse  problem.  Hence,  it  is  natural  to  solve  it  via  the  minimization  of  an  objective 
functional  J[q,  f]  with  respect  to  q.  In  our  case  the  functional  J[q,  f]  consists  of  two  functionals, 

J  =  J f  +  Js-,  (3-1) 

where  Jf  estimates  the  fulfillment  of  the  conditions  (1.6)  and  (1.7)  and  Js  estimates  fulfillment  of  condition 
sup  q  £  D. 

The  functional  Jf  consists  of  four  functionals 

Jf  =  Jh\  +  Jgi  +  Jh2  +  Jg2i  (3-2) 

where  functionals  J/M ,  Jgi,  Jh2 ,  Jg2  respectively  estimate  fulfillments  of  conditions  u|r1T  =  hi,  ux |r1T  =  g  1, 
u\r2T  h‘2 ,  uy\r2T  g2 • 
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The  functional  Js  has  the  form 


Js[q}  =  (Ml~h\\2D)/2  =  (q,Ssq)/2, 
where  Ss  is  a  symmetric  non-negative  matrix  with  elements 

Ss,nmik  =  (a  +  T)(b  +  T)8nlSmk  -  X^k. 
Introduce  operators  A /ll ,  Agi ,  Ah2 ,  Ag2  such  that 

Ahlq  =  (Aq)  |  Agiq=(Aq)x | 


(3-3) 


(3-4) 


Ah2q  =  (Aq) 


r2T’ 


Ag2  q  =  (Aq) 


(3.5) 


y  ir2T  * 


Let  ||  •  ||  denotes  the  Z/2-norm.  We  can  write  the  expression  for  J/ll  as  follows 


JhAQJ]  =  P/uQ-  M2/2  =  PftiQll2/2  -  (Ahlq,  hi)  +  \\hi\\2/2  =  ( q,Shlq)/2  -  (q,hi)  +  ||/ii||2/2,  (3.6) 

where  5/ll  =  A*hiAkl  is  the  symmetric  non-negative  matrix;  h\  =  A*hJii  is  the  vector,  and  (v,w)  denotes 
the  scalar  product, 

( V,W )  —  ^  '  ^nmU'nm • 
nm 

The  expressions  for  elements  of  the  matrix  S), ,  and  the  vector  h\  depend  on  the  problem  we  consider, 
i.e.,  the  (^-problem  or  the  ■(/(-problem.  We  mark  the  elements  for  the  (^-problem  by  the  superscript  and 
the  elements  for  the  '(/’-problem  by  the  superscript  ■(/>.  Then 


c + 

°h\,nmlk 


=  SnSlY^^mlk, 


a b 

°hi,nmlk 


rb+T  n  i 

Knm  =  sn  /  Ym(y )  /  cos  (u)nmt)  hi(y ,  t)  dt  dy, 

Jo  Jo 


foSS 

—  Q  ^nmlk 

—  ■ 

UnmUlk 


+  = 

ul,nm 


rb+T 


OJnm  Jo 


Ym(y )  /  s\n(unrnt)hi(y,t)dtdy , 


where 


sn  =  Xn(0)  =  sin  (7rn/2); 

[ZYm(y)Yk(y)dy- 

Jo 

®nmlk  =  /  cos  pnmt)  COS  (Wifct)  dt; 

Jo 

$nmlk  =  /  sin  (Vnmt)  sin  (uikt)  dt. 

Jo 

We  obtain  analogous  expressions  for  other  functionals.  We  get  for  Jgi 

Jgihf]  =  \\AgiQ  -  dill/2  =  ( Q,Sgiq)/2  -  ( q,gi )  +  ||5i|j2/2, 


where  Sgi  =  A*giAgi  is  the  matrix  with  elements 


a f 

gi,nmlk 


=  cnqP+T$cc 


ir2nl 


mk  nmlk  _|_  j^2  > 


cb  =  „  „,W>+T 
°gi,nmlk  i'"n^1mk 


^nm^lk  4((2  + 


(3.7) 

(3.8) 

(3.9) 
(3.10) 

(3-11) 

(3.12) 

(3.13) 
(3-14) 

(3.15) 

(3.16) 

(3.17) 
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(3.18) 


cn  =  -X^(O)  =  cos  (7rn/2); 


gi  =  A*t  g\  is  the  vector  with  components 


ircnn  rb+T  rT 

9i,nm  =  2(q  "  T )  Ym{y)  cos  (unmt)gi(y,  t)  dt  dy, 


-0  _  7T Cn7l 

,nm 


rb+T 


2 (cl  +  T)cjnrn  J q 


Ym(y)  /  sin  (a;  nm  t)gi(y,t)dtdy. 


We  obtain  for  J^2 


<42[<?]  =  -  /i2||2/2  =  (■ q,Sh2q)/2  -  ( q,h2 )  +  \\h2\\2 /2, 


where  entries  of  the  matrix  Sh  =  are 


qT  —  q  q,  XaJrT <bcc 
°h2,nmlk  ~  bm^k^ni  ^ nmlki 

<f,SS 

Qip  _  „  vd+T  ^ nmlk 

°h2,nmlk  ~  bmbk^ni 


^nm^lk 


ra+T 

Xni=  /  Xn(x)Xi(x)  dx; 

Jo 

h2  =  A*hoh2  is  the  vector  with  components 

T 

Xn(x)  /  cos  (. (jjnmt)h2{x ,  t)  dt  dx, 

Jo 


h2,nm  —  s m 


ra+T 


We  obtain  for  J, 


92 


ra+T  rT 

l2nm  =  sm  /  Xn(x)  /  sin  {u>nmt)h2(y,  t)  dt  dy. 

Jo  Jo 


J92+  /]  =  P32Q  -  52|!2/2  =  (Q,  S92q)/2  -  ( q,g2 )  +  \\g2\\2/2, 


where  Sg2  =  A*0Ag2  is  the  matrix  with  entries 


Sl,nmik  = 


nl  ^nmlk 


ir2mk 
4(6 +  T)2’ 


ir2mk 


a+T  nmlk 


^ Ho  '/?  77?  /  A*  C?7?  C/(‘  JJ  1  ,  ^  .  - 

92’  UnmUlk  4(6 +  T)2 


g2  =  A*2g2  is  the  vector  with  components 


t  cmm 

92, nm  — 


ra+T 


2(6 +  T)  y0 


-X'n(x)  /  cos  ((jJnmt)g2(y,  t)  dt  dy, 


/  0 


xfj  _  7T  cmra 

52-nm  “  2(6  + tv 


/*a+T  rT 

/  Xn(x)  /  sin  (a; 

nm  t)g2{y,t)  dt  dy. 

nm  Jo  Jo 


Using  (3.6),  (3.15),  (3.21),  (3.27),  we  obtain  the  following  expression  for  J/: 

Jf  =  prq  -  /f/2  =  (q,  S,q)/2  -  (q,  /)  +  ||/||2/2, 

where  Ar  is  the  linear  operator:  Arq  =  {Ahlq,  Agiq,  Ah2q,  Agoq},  f  =  {hi,  gi,  h2,  g2}  i 
Sf  =  A*rAr  is  a  symmetric  non-negative  matrix 


sf  =  shl  +  sgi  +  sh2  +  s. 


92  1 


(3.19) 

(3.20) 

(3.21) 

(3.22) 

(3.23) 

(3.24) 

(3.25) 

(3.26) 

(3.27) 

(3.28) 

(3.29) 

(3.30) 

(3.31) 

(3.32) 

a  vector  function. 

(3.33) 
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and  /  =  A^f  denotes  the  vector 


f  =  h1  +  gi  +  h2  +  g2- 


(3.34) 


Combining  (3.32)  and  (3.3),  we  obtain 

J[<h  f]  =  JfW\  +  Js[q\  =  (9,  Sq)/2  -  ( q ,  /)  +  H/ll/2,  (3.35) 

where  S  is  a  symmetric  non-negative  matrix  with  elements 

S  =  Sf  +  Ss.  (3.36) 

It  follows  from  (3.35)  that  J  is  a  quadratic  functional.  It  is  well  known  that  one  of  the  most  effective 
methods  of  minimization  of  quadratic  functionals  is  the  conjugate  gradient  method.  Thus,  we  will  use  it. 
This  method  can  be  schematically  described  as  follows. 

1.  On  the  fc-th  iteration  we  compute  the  gradient  of  J  with  respect  to  q 

9k  =  VqJfc  =  Sqk  -  f.  (3.37) 

Suppose  that  \gk\  <  e,  where  e  is  an  a  priori  chosen  small  number.  Then  we  stop  the  iterative 
process  and  consider  the  vector  qk  as  the  solution.  Otherwise,  we  go  to  the  next  step. 

2.  Compute  the  conjugate  gradient 

Pk  =  9k  -  PkPk-i,  Pk  =  (■ Spk-i,gk)/(Spk-i,Pk-i )•  (3.38) 

For  k  =  0  set  po  =  g o- 

3.  Calculate  the  step  of  descent 

ak  =  {Pki9k)/{Spk,Pk)-  (3.39) 

4.  Compute  the  new  vector 

9k+ 1  =  9k~  akpk  (3.40) 

and  return  to  the  first  step. 


It  well  known  that  the  number  of  iterations,  which  is  required  to  reach  the  minimum  of  the  quadratic 
functional  equals  to  the  dimension  of  the  vector  q  [14] . 

To  estimate  the  accuracy  of  our  solution,  we  refer  to  the  Tikhonov  concept  of  solutions  of  ill-posed 
problems,  see  Tikhonov  and  Arsenin  [13].  By  this  concept  one  needs  to  assume  the  existence  of  an  ” ideal” 
exact  solution  with  an  ideal  exact  data.  Then  one  should  assume  that  a  certain  error  is  perturbing  the 
ideal  data,  which  leads  to  a  non-ideal  data,  and  estimate  the  difference  between  the  computed  solution 
corresponding  that  non-ideal  data  and  the  exact  solution.  Thus,  we  arrive  at  the  following  lemma 

Lemma  1.  Let  q*  be  the  exact  solution  of  the  inverse  problem  with  the  exact  data  /*  =  (L*,q*)  and  q 
be  the  minimizer  of  the  functional  J\q,  f]  with  the  data  f  =  ( h,g ).  Then 

(^/, max) 1/2  || /  f*\\  +  II i  /]  II 

^min 

Wf-Ml  +  V^AoJ] 

\J  ^min 

where  Sf>max  is  the  maximal  eigenvalue  of  the  matrix  Sf  and  smjn  is  the  minimal  eigenvalue  of  the  matrix  S. 


(3.41) 


Proof.  Since  g*  is  the  exact  solution  with  the  exact  data  we  have  J[q*,f*]  =  0  and  Vq  J[q*,  fast]  = 
Sq*  —  f*  =  0.  Therefore 

V,J[g,  /]  =  (Sq  -  f)  -  ( Sq *  -  /*)  =  S(q  -  q*)  -  (/  -  /*). 

Using  the  triangle  inequality,  we  obtain 

||S(g-<Z*)ll  <  ||/-/*ll  +  l|VqJ[g,/]||. 

Taking  into  account  the  inequalities 

II  f  ~  f*  II  =  Pf  (/  -  /*)||  <  V^/.m ax  11/  -  /*ll> 

PCq  -  q*)||  >  Smin||q  -  <7*  || , 

we  obtain  (3.41). 

To  prove  (3.42),  consider  the  gradient  of  J  with  respect  to  /, 

VfJ[q,f\  =  f-Arq. 

Since  /*  =  A^q*,  we  have 

V/J[g,/]  =  (/  -  P<?)  -  (/*  -Avq*)  =  (/  -  /*)  -  ^r(gi  -  <?2)- 
Using  the  triangle  inequality,  we  obtain 


Pr(g  -  <?*)H  <  ||/  -  /*||  +  ||V/J[g,/]||, 

which  with  relations 

\\Ar(q  -  q*)  ||  >  Vsmin  ||<7  -  g*||, 

\\VfJ[q]\\  =  \\f-Arq\\  =  y/2J\q,f], 

yields  (3.42). 


Remark  1.  Note  that  the  gradient  of  the  objective  functional  with  respect  to  the  data  is  used  in  the 
proof  of  this  lemma.  The  authors  are  not  aware  about  other  publications  which  would  employ  this  idea. 


4  Numerical  results 


We  solve  numerically  both  ip-  and  ^-problems  for  the  test  function  q*(x,y)  given  by 


q*(x,y) 


sin  (2irx/a)  sin  (2ny/b),  (x,y)  £  D; 

0,  (x,y)#D. 


The  following  values  of  parameters  are  taken:  a  =  1,  b  =  1,  T  =  1.5,  i.e.  D  =  (0,1)  x  (0,1),  U  = 
(—2.5,  2.5)  x  (—2.5, 2.5)  (see  Figure  2).  The  domain  U  is  partitioned  in  M  =  26  =  64  equal  parts  with 
respect  to  the  variable  y  and  in  N  =  26  =  64  equal  parts  with  respect  to  x.  The  time  interval  (0,  T)  is 
partitioned  in  K  =  100  equal  parts. 

First,  we  solve  the  direct  problem  and  compute  the  functions  {/ii*,  gu,  h2*,g2*}  =  /*•  Figure  3  illustrates 
the  functions  /12*,  <?2*  (^1*  =  9i*  =  92*  due  to  the  symmetry)  for  ip-  and  ^-problems  respectively. 

Next,  we  solve  the  inverse  problem  with  the  data  /  =  /*.  Each  of  the  integrals  entering  the  expressions 
for  entries  of  the  matrices  S^,  Sgi,  Sh.2,  Sg2 ,  Ss  and  vectors  hi,  g\,  h-2,  92  (see  the  previous  section)  is 
approximated  by  the  sum 

f  u(x)  dx  &  Ax  u(n  Ax), 
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Figure  2:  The  test  function  q*(x,y) 


Figure  3:  Solution  of  the  direct  problem:  h,2*(x,t)  (left);  g2*(x,t)  (right) 


where  Ax  is  the  discretization  step.  Using  q  =  0  as  the  initial  guess,  we  proceed  with  the  conjugate  gradient 
method  as  long  as  either  the  number  of  iterations  is  less  than  300  or  the  norm  of  the  gradient  of  the  functional 
is  exceeds  10~10. 

The  results  of  the  numerical  solution  of  <p-  and  ^-problem  are  presented  on  Figures  4,  5  respectively. 
One  can  see  that  in  order  to  solve  the  ^-problem,  only  30  iterations  are  required  to  reach  ||V9J||  <  10-10. 
Whereas  the  norm  of  the  gradient  of  the  functional  for  the  ^-problem  ||VJg||  ~  2  •  10~2  after  300  iterations. 
To  explain  this  difference,  we  consider  the  eigenvalues  of  matrices.  The  results  of  numerical  computations 
of  the  eigenvalues  are  presented  in  Table  1.  It  follows  that  the  condition  number  (x  =  smax/.smin)  of  the 
matrix  Sf  +  Ss  for  the  (^-problem  is  yf  =  8903/6.536  ~  1362,  while  the  condition  number  for  the  problem 
is  =  9.605/0.832  ~  11.5.  Therefore,  yf  almost  100  times  greater  than  y^ ,  which  explains  why  the 
convergence  rate  for  the  ^/-problem  is  much  better  than  for  the  (^-problem. 


Remark  2.  As  it  was  mentioned  above  in  the  description  of  the  algorithm  for  the  inverse  problem,  to 
achieve  the  minimum  it  is  required  dim  q  iterations.  In  our  case  dim  q  =  NM  =  64  •  64  =  4096.  However,  we 
achieve  the  minimum  after  about  100  iterations,  which  is  a  very  good  number.  We  cannot  yet  fully  explain 
this. 
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Figure  4:  Solution  of  the  (^-problem  (left)  and  dependencies  of  the  norm  of  the  gradient  of  func¬ 
tional  ||V?J[^,/*]||  (solid  curve),  the  functional  J[ipki  f*\  (dashed  curve),  the  norm  of  deviation  from  the 
exact  solution  || —  <p*||  (dash-dotted  curve)  on  the  iteration  number  k  (right) 


Figure  5:  Solution  of  the  ^-problem  (left)  and  dependencies  of  the  norm  of  the  gradient  of  func¬ 
tional  ||Vq J[if)k,  /*] ||  (solid  curve),  the  functional  J[^fc,/*]  (dashed  curve),  the  norm  of  deviation  from  the 
exact  solution  \\ipk  —  t/>*||  (dash-dotted  curve)  on  the  iteration  number  k  (right) 


Table  1:  The  eigenvalues  of  matrices 


s 

^min 

^max 

Ss 

0 

6.25 

S 7 

1.861  •  10~4 

8900 

sf 

7.988- 10-7 

3.455 

sf  +  ss 

6.536 

8903 

S'f  +  Ss 

0.832 

9.605 
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Table  2:  The  eigenvalues  and  condition 
number  of  the  matrix  Sv  (7) 


7 

^min 

°  max 

X 

0 

1.861  •  10"4 

8900 

4.78  •  107 

1 

6.536 

8903 

1360 

5 

16.034 

8914 

556 

10 

21.104 

8928 

423 

30 

27.569 

8984 

326 

50 

29.782 

9041 

304 

100 

32.032 

9191 

287 

200 

33.567 

9518 

284 

300 

34.206 

9886 

289 

1000 

35.320 

13563 

384 

10000 

35.881 

69547 

1940 

100000 

35.945 

628557 

17500 

Table  3:  The  eigenvalues  and  condition 


number  of  the  matrix  S ^  (7) 


7 

^min 

^max 

X 

0 

7.988 -lO"7 

3.455 

4.325  •  107 

0.1 

0.281 

4.023 

14.27 

0.2 

0.498 

4.638 

9.304 

0.3 

0.628 

5.257 

8.360 

0.4 

0.692 

5.879 

8.496 

0.5 

0.734 

6.501 

8.853 

0.6 

0.764 

7.123 

9.312 

0.7 

0.787 

7.745 

9.828 

1 

0.832 

9.605 

11.53 

5 

0.931 

34.24 

36.75 

10 

0.946 

65.77 

69.52 

100 

0.959 

627.8 

654.4 

Figure  6:  Solution  of  the  (^-problem  with  7  =  200  (left)  and  dependencies  of  the  norm  of  the  gradient  of 
functional  ||  Vq  J[ipki  /*]  II  (solid  curve),  the  functional  /*]  (dashed  curve),  the  norm  of  deviation  from 
the  exact  solution  \\ijjk  —  i/bll  (dash-dotted  curve)  on  the  iteration  number  k  (right) 


2.5  K 


Figure  7:  Solution  of  the  •i/’-problem  with  7  =  0.3  (left)  and  dependencies  of  the  norm  of  the  gradient  of 
functional  || Vq  Jfyk,  /*]  ||  (solid  curve),  the  functional  J[ipk,  f*]  (dashed  curve),  the  norm  of  deviation  from 
the  exact  solution  W'tpk  —  i/>*||  (dash-dotted  curve)  on  the  iteration  number  k  (right) 
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To  decrease  the  condition  number  x,  we  consider  the  following  functional 


J(7)  =  J/  +  7JS,  (4.1) 

where  7  is  a  parameter.  Dependencies  of  eigenvalues  and  the  condition  number  of  the  matrix  5(7)  =  Sf+'ySs 
on  the  parameter  7  for  cp-  and  (/-problems  are  presented  in  Tables  2  and  3  respectively.  One  can  see  that 
the  minimal  condition  number  Xmm  =  284  for  the  (^-problem  is  reached  when  7%n  =  200,  while  for  the 
'(/’-problem  Xmin  =  8-36  is  reached  when  7^in  =  0.3. 

Figure  6  illustrates  the  numerical  solution  of  the  ^-problem  with  7  =  200,  and  Figure  7  illustrates 
the  numerical  solution  of  the  (/-problem  with  7  =  0.3.  The  convergence  rate  of  the  (/9-problenr  is  notably 
improved,  since  only  106  iterations  is  required  now  to  reach  ||VqJ[g,  f]\\  <  10~10.  The  convergence  rate  of 
the  ■0-problem  did  not  change,  because  the  condition  number  for  7  =  1  does  not  differ  significantly  from  the 
condition  number  for  7  =  0.3  (see  Table  3). 

5  The  case  of  a  random  noise  in  the  data 

To  estimate  the  influence  of  the  random  noise  in  the  data  /  on  the  solution  of  the  inverse  problem,  we 
introduce  this  noise  in  the  data  and  solve  the  resulting  inverse  problem.  So,  the  precise  data  /*  is  replaced 
with  the  noisy  data  f£  as 

fe,kn  =  (1  +  ernd  ())/*!fcn-, 

where  e  is  the  noise  level  and  rnd  ()  is  the  result  of  the  random  number  generator  with  the  normal  distribution 
in  [-1,1]. 

Tables  4,  5  illustrate  the  results  of  numerical  solutions  of  both  ip-  and  (/- problems  for  various  levels  e  of 
the  noise.  Since  ||<^*||  =  ||(/*||  =  1/2,  then  the  relative  differences  can  be  obtained  from  norms  || ip  —  <^*||, 

1 1  (/  —  (/*]  |  given  in  these  tables  via  the  multiplication  by  5.  In  addition  to  the  norm  of  the  deviation  of  the 
computed  solution  from  the  exact  one  ||q  —  q*  ||,  we  present  the  norm  of  the  deviation  of  the  noisy  data  from 
the  exact  one  \\f£  —  /*||,  the  value  of  the  functional  J[q\  and  estimates  (3.41)  and  (3.42).  The  norm  of  the 
gradient  ||V9J||  <  10~10  at  the  last  iteration  for  all  e  and,  therefore  is  not  presented.  The  values  of  s/)max> 
sm;n  entering  estimates  (3.41),  (3.42)  are  taken  from  Table  1  and  Tables  2,  3  {pf  =  200,  7^  =  0.3).  An 
interesting  feature  of  these  results  is  that  even  at  e  =  0.5,  which  corresponds  to  the  50%  random  noise  in  the 
data  the  relative  difference  between  computed  and  exact  solutions  does  not  exceed  4%  for  the  (^—problem 
and  10%  for  the  (/—  problem. 

One  can  see  that  the  estimate  (3.42)  is  more  accurate  than  (3.41).  To  explain  this,  consider  the  case  when 
the  number  of  iterations  is  not  limited.  In  this  case  we  would  obtain  such  a  solution  q  that  ||VgJ[g, f£]  ||  =  0 


Table  4:  Deviation  of  numerical  solution  for  99-problenr  with  =  sin  (2irx/a)  sin  (27T y/b) 


£ 

II A  -  All 

J  [<P\ 

\w 

-  ¥>*ll 

Estimate  (3.41) 

Estimate  (3.42) 

0 

0 

3.55-  10~15 

1.24 

10-13 

2.38  ■  10“12 

1.45  •  10-8 

0.001 

0.001 

8.27-  10-7 

1.29 

■  10-5 

0.003 

4.55-  10"4 

0.005 

0.006 

2.09-  10-5 

6.35 

■  10-5 

0.018 

0.002 

0.01 

0.013 

8.39-  10-5 

1.88 

■  10“4 

0.038 

0.004 

0.05 

0.066 

2.01  ■  10~3 

6.87 

■  10“4 

0.187 

0.022 

0.1 

0.134 

8.17-  10~3 

1.53 

■  10-3 

0.377 

0.045 

0.5 

0.682 

0.210 

7.99 

■  10-3 

1.918 

0.229 

Table  5:  Deviation  of  numerical  solution  for  ■(/’-problem  with  (/*  =  sin  (27rx/a)  sin  (2ny/b) 


£ 

ll/e  -  All 

J[ip\ 

IIV’  - 

V>*|l 

Estimate  (3.41) 

Estimate  (3.42) 

0 

0 

5.551 

10~17 

3.865  • 

10”11 

7.45-  10"11 

1.32  -  10~8 

0.001 

1.57  ■  10“4 

1.129 

•  10“8 

3.408  ■ 

10“5 

4.64-  10~4 

3.87-  10~4 

0.005 

7.56  ■  10"4 

2.595 

•10-7 

1.672- 

10“4 

2.23-  10~3 

1.86  •  10'3 

0.01 

1.51  ■  10-3 

1.047 

•  10“6 

3.215- 

10“4 

4.48-  10-3 

3.73  •  10~3 

0.05 

7.76  ■  10“3 

2.709 

■  10“s 

1.854  - 

10“3 

0.022 

0.019 

0.1 

0.015 

1.111 

■  10“4 

3.136- 

10-3 

0.046 

0.038 

0.5 

0.076 

2.591 

•  10“3 

0.019 

0.225 

0.187 
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Figure  8:  Solution  of  the  (^-problem  with  e  =  0.5  (left)  and  dependencies  of  the  norm  of  the  gradient  of 
functional  ||  Vq  J[ipk,  fe]  ||  (solid  curve),  the  functional  J[<pk,  fe]  (dashed  curve),  the  norm  of  deviation  from 
the  exact  solution  || ip^  —  p* ||  (dash-dotted  curve)  on  the  iteration  number  k  (right) 

and  J[q,  fe]  «  ||  fe  -  f  ||2/2.  Then  _ 

^/,max 
^min 

Therefore,  if  s/,max/smin  >  4,  then  the  estimate  (3.42)  is  more  accurate  than  (3.41).  However,  if  s/,max/smin  < 
4,  then  the  estimate  (3.41)  is  more  accurate  than  (3.42).  Since  s/,max/smin  >  4  for  both  ip-  and  i/j-problems, 
then  the  estimate  (3.42)  should  give  lesser  values  than  (3.41).  It  follows  from  Tables  4,  5  that  the  computed 
deviation  of  the  solution  is  almost  10  times  less  than  one  predicted  by  the  estimate  (3.42).  This  may  be 
explained  by  the  smoothness  of  the  exact  solution  y* . 

We  now  consider  the  (^-problem  for  the  (5-like  test  function  p*{x,  y)  =  5(x— 0.9)S(y— 0.9).  This  is  an  exact 
analog  of  the  problem  of  time  reversal.  The  discrete  analog  of  S(x— 0.9)6(y— 0.9)  is  <p*tnm  =  ^nn'^mm' / y/hxhy, 
where  5nn/  =  0  if  n  =  n'  and  6nni  =  0  otherwise;  n'  is  the  index  of  the  grid  point  close  to  x  =  0.9;  m!  is 
the  index  of  the  grid  point  close  to  y  =  0.9;  hx  and  hy  are  discretization  steps  along  the  x-  and  y-axis 
respectively.  Figure  8  illustrates  the  solution  of  the  problem  with  the  quite  high  noise  level  e  =  0.5.  The 
norm  of  the  deviation  || p  —  p*\\  =  2.78  •  10-2,  while  the  estimate  (3.42)  gives  2.31  •  10-1.  Hence,  we  have 
now  ~  7%  of  the  relative  error  in  the  solution  at  50%  of  the  random  noise  in  the  data.  Thus,  the  high 
accuracy  of  numerical  solutions  cannot  be  attributed  only  to  the  smoothness  of  the  exact  solutions.  So,  we 
provide  some  insights  of  this  phenomenon  in  the  next  section.  We  note  that  Figure  8  demonstrates  a  quite 
good  refocusing  of  the  time  reversed  wave  field  in  the  presence  of  a  large  amount  of  random  noise  in  the 
data. 

6  Summary 

We  are  motivated  by  an  interesting  experimentally  observed  phenomenon  of  refocusing  of  time  reversed 
wave  fields  [3].  Using  a  mathematical  model  of  time  reversal  proposed  in  [6],  we  have  developed  method 
for  the  numerical  solution  of  the  inverse  problem  of  determining  an  unknown  initial  condition  in  the  2-D 
wave  equation  utt  =  uxx  +  uyy .  The  support  of  the  unknown  initial  condition  is  in  the  first  quadrant  of  the 
(x,  y)—  plane  and  the  lateral  Cauchy  data  are  given  at  finite  parts  of  coordinate  axis.  A  number  of  numerical 
results  is  presented,  including  the  one  (Figure  8),  which  directly  models  the  phenomenon  of  refocusing  in  the 
time  reversal  experiment.  We  have  estimated  convergence  rates  and  accuracy  of  solutions  of  our  algorithm 
for  both  <p -  and  ^—problems  both  analytically  and  numerically.  An  interesting  feature  of  our  computational 
results  is  their  surprisingly  high  accuracy  even  in  the  presence  of  a  large  amount  of  random  noise  in  the 
data.  Indeed,  at  50%  noise  in  the  data  the  relative  error  of  our  solutions  repeatedly  did  not  exceed  a  few 


q-  q*  1 1  (3.4i)  _  l  / 
Q  —  Q*  II  (3.42)  2  V 
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Figure  9:  Solution  of  the  (^-problem  with  e  =  0.3  (left)  and  e  =  0.5  (right) 

per  cent.  In  our  opinion,  this  can  be  explained  by  the  existence  of  a  priori  Lipschitz  stability  estimate  [6] 
for  this  inverse  problem.  We  recall  that  the  Lipschitz  stability  estimate  is  the  best  possible.  Observe  that 
a  similar  high  accuracy  in  the  presence  of  a  large  amount  of  noise  in  the  data  was  also  demonstrated  in 
the  much  earlier  work  [10],  where  the  2-D  wave  equation  was  solved  with  the  lateral  Cauchy  data  at  the 
boundary  of  a  square.  Because  this  square  is  a  bounded  domain,  then  the  Lipschitz  stability  of  the  problem 
of  [10]  takes  place  [9,  7,  5],  which  is  similar  with  the  case  of  the  above  Inverse  Problem  2. 

Results  of  this  paper  provide  a  certain  computational  confirmation  of  the  high  degree  of  refocusing  of 
time  reversed  wave  fields  in  the  presence  of  “ reflecting  boundaries  as  waveguides  or  reverberating  cavities” 
(p.  R2  in  [3]),  while  Lipschitz  stability  estimates  of  [6,  9,  7,  5]  might  be  viewed  as  ones  providing  a  theoretical 
foundation  of  this  experimental  observation.  On  the  other  hand,  Figure  9  displays  solution  of  the  ^-problem 
for  yj*(x,y)  =  6(x  —  0.9 )5(y  —  0.9)  (i.e.,  for  the  same  99*  as  above)  and  for  the  case  when  the  Cauchy  data 
are  given  only  at  the  y-axis,  i.  e.  at  the  line  Fi  (see  Figure  1).  We  have  taken  e  =  0.3  and  0.5,  which  is  30% 
and  50%  the  random  noise  in  the  data.  We  got  \\<p  —  ip*\\  =  0.36  and  1.02,  or  90%  and  250%  of  relative 
error.  It  is  clear  that  this  result  is  much  worse  than  one  of  Figure  8,  although  the  noise  level  is  less  than  in 
the  previous  case.  Although  an  analog  of  the  stability  estimate  of  [6]  (i.e.,  in  the  “continuous  case”)  is  not 
yet  proven  for  the  geometry  corresponding  to  Figure  9,  but  numerical  tests  of  Figure  9  indicate  that  such 
an  estimate  is  likely  much  weaker  than  the  Lipschitz  estimate,  because  of  different  geometries 
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